Lyapunov exponents from unstable periodic orbits in the FPU-/? model 
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In the framework of a recently developed theory for Hamiltonian chaos, which makes use of 
the formulation of Newtonian dynamics in terms of Riemannian differential geometry, we obtained 
analytic values of the largest Lyapunov exponent for the Fermi-Pasta-Ulam-/3 model (FPU-/3) by 
computing the time averages of the metric tensor curvature and of its fluctuations along analytically 
known unstable periodic orbits (UPOs). The agreement between our results and the Lyapunov 
exponents obtained by means of standard numerical simulations supports the fact that UPOs are 
reliable probes of a general dynamical property as chaotic instability. 
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INTRODUCTION 



Unstable periodic orbits (UPOs) are widely studied in 
the framework of classical nonlinear dynamical systems 
, since they form the " skeleton" Q of the phase space 
of these systems and are very sensitive to local character- 
istic features of the dynamics. Among the applications 
of the studies of UPOs we can mention: a characteriza- 
tion of dynamical systems 13, control of classical chaos, 
semiclassical quantization [J. Furthermore, they are use- 
ful for a characterization of quantum chaos and for the 
description of some thermodynamical properties of dy- 
namical systems with many degrees of freedom. 

The present paper aims at lending further credit to 
the common wisdom of relevance of UPOs for chaotic 
dynamics. Our contribution to this subject stems from a 
Riemannian geometric approach to the study of Hamil- 
tonian chaos. 

It is well known that the degree of chaoticity of a dy- 
namical system is measured by the largest Lyapunov 
exponent Ai which provides an average dynamical in- 
stability growth rate in terms of the local growth rate 
of the distance of nearby trajectories, averaged along a 
sufBciently long reference trajectory. The largest Lya- 
punov exponent Ai for standard Hamiltonian systems, 
described by Hamiltonian functions of the form H = 
SiLi \Pi + ^(91 7 • • • ' 9^)7 is computed by numerically 



integrating the tangent dynamics equation 
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along a reference trajectory q{t) = [qi{t), ..,qisi(t)], 
and then Ai = Yant^^l/2t\Qg{Y.f^^[if{t) + 
/Sill [4^(0) + Cl(0)]). In the conventional theory 
of chaos, dynamical instability is caused by homoclinic 
intersections of perturbed separatrices, but this theory 
seems not adequate to treat chaos in Hamiltonian 
systems with many degrees of freedom. In this case, the 
direct numerical simulation is the only way to compute 
Ai. 

Recently, it has been proposed by Pettini j5| to tackle 
Hamiltonian chaos in a different theoretical framework 
with respect to that of homoclinic intersections. This 
new method resorts to a well known formulation of 
Hamiltonian dynamics in the language of Riemannian 
differential geometry: the mechanical trajectories of a 
dynamical system can be viewed as geodesies of a Rie- 
mannian manifold endowed with a suitable metric. In 
this framework, it is possible to relate the instability 
of a geodesies flow with the curvature properties of the 
underlying "mechanical" manifold through two geomet- 
ric quantities: the Ricci curvature and its fluctuations. 
These two geometric quantities, in principle averaged 
along a generic geodesic, enter a formula, derived by Pet- 
tini et al. in Ref.Q, which allows the analytic compu- 
tation of the largest Lyapunov exponent for a generic 
Hamiltonian system. However, since the mentioned time 
averages are in general not analytically knowable, one 
has to replace them with microcanonical averages which 
coincide with time averages when the number of degrees 
of freedom is large and the dynamics is chaotic. In fact. 
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under these circumstances, the measure of regular orbits 
in phase space - at physically meaningful energies - is 
vanishingly small, thus the dynamics is bona fide ergodic 
and mixing. 

Therefore, the analytic computation of the largest Lya- 
punov exponent can be done whenever the simplifying 
hypotheses of Ref.|^ are justified and the microcanonical 
averages of the mentioned geometric quantities are ana- 
lytically computable. This is just the case of the FPU-/3 
model which has been considered in Ref. 0|. 

It is the purpose of the present work to show that, in 
some special case, the above mentioned replacement of 
time averages with microcanonical ones can be avoided: 
provided that the time averages of the Ricci curvature 
and of its fluctuations are analytically computed along 
some unstable periodic orbits, a reasonable analytic es- 
timate of the values of Ai can be obtained without re- 
sorting to microcanonical averages. It is somewhat sur- 
prising, and undoubtely very interesting, that unstable 
periodic orbits make something like an "importance sam- 
pling" of the relevant geometric features of configuration 
space which are needed to estimate the average degree of 
chaoticity of the dynamics, measured by Ai. A similar 
problem was already addressed in 0, where the authors 
gave an analytical estimate of the largest Lyapunov expo- 
nent at high energy density for the Fermi-Pasta-Ulam-/3 
model by computing the average of the modulational in- 
stability growth rates associated to unstable modes. 



2. GEOMETRY AND DYNAMICS 

Let us summarize the geometrization of Newtonian dy- 
namics tackled in . It applies to standard autonomous 
systems described by the Lagrangian function (all the 
indices run from 1 to degrees of freedom) 

^(9''j) = ^E«^'^-('?M''?'-^(9) , (2) 

ik 

where a^fc is the kinetic energy tensor that in terms of 
the total energy E and kinetic energy, reads 

J2a^kq'q''^2iE-V)^2W , (3) 

ik 

Following the method due to Eisenhart jSj, the differen- 
tiable A^— dimensional configuration space A4, on which 
the lagrangian coordinates (g^,...,g^) can be used as 
local coordinates, is enlarged. The ambient space thus 
introduced embodies the time coordinate and is given 
as X R^, with local coordinates g^, . . . , g^, q^~^^), 
where {q^, . . . ,q^) € M, q° S M is the time coordinate, 
and g R is a coordinate closely related to Hamilton 

action. With Eisenhart we define a pseudo-Riemannian 



non-degenerate metric g^^ on A4 x as 
dsl ^ fM-^ dqf" (g) dq"" = dq° dq^+^ + dq^+^ ® dq° + 
Qij dq"- (g) dq' - 2V{q) dq° (K) dq° . 

ij 

(4) 

Natural motions are now given by the canonical projec- 
tion TT of the geodesies of {A4 x M^,^^;) on the config- 
uration space-time: tt : x ^ x M. However, 
among all the geodesies of g,^ the natural motions belong 
to the subset of those geodesies along which the arclength 
is positive definite 

ds^ = g^.^dq^dq" = 2C'^df > , (5) 

flu 

where C is a real arbitrary constant. More details can 
be found in 

The stability of a geodesic flow is studied by means 
of the Jacobi— Levi-Civita (JLC) equation for geodesic 
spread. In local coordinates and in terms of proper time 
s the JLC equation reads as 

ijr 

where J is the Jacobi vector field of geodesic separation, 
where the covariant derivative is given by VJ'^/ds = 
dJ'^/ds + ^ij^ij dq^ /dsJ^ , and R^j^ are the compo- 
nents of the Riemann-Christoffel curvature tensor which, 
in terms of the Christoffel coefficients Fjfj, are 

t 

where dj = d/dq^ . The Christoffel coefficients, in turn, 
are defined as 

^jk = ^ XI 5 ™ idj9km + dkgmj - drag^k) ■ (8) 

m 

In 1^ it has been shown that the Jacobi equation 10, 
written for the Eisenhart metric of the enlarged configu- 
ration space, nicely yields the standard tangent dynam- 
ics equation Moreover, under suitable simplifying 
hypotheses, mainly of geometric type, in Ref. [6| it has 
been shown that Eq.® can be replaced by a scalar ef- 
fective equation 

^ + (fcfl),^ + -^l={s'Kj^yjM^)^ = , (9) 

where ip stands for any of the components J* of the Jacobi 
field, since in this effective picture all of them obey the 
same equation. Moreover, Kr = J2ijk 9^'' ^\kj 
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Ricci curvature and ku = Kn/^N — 1) which, for the 
Eisenhart metric, takes the simple form 



kR{q) 



AV 



1 " 
N ^ 



(10) 



In Eq.lPJ, rj{s) is a gaussian white noise with zero mean 
and unit variance, and {■)s stands for time averaging 
along a reference geodesic. Time averages {kfi)s and 
{6'^Kji)s of Ricci curvature and of its second moment re- 
spectively, cannot be known analytically for a chaotic or- 
bit, hence the need of an assumption of ergodicity allow- 
ing the replacement of time averages by microcanonical 
averages on a constant energy surface T,e, correspond- 
ing to the energy value E of interest. At variance with 
time averages along chaotic orbits, microcanonical aver- 
ages can be computed analytically for some models. It is 
worth remarking that, after the replacement of time aver- 
gaes by means of static microcanonical averages (fc^)^^ 
and ((J^iffl)^^, the scalar equation @ is independent of 
the numerical knowledge of the dynamics. 

Then the largest Lyapunov exponent for the effective 
model given by Eq. defined as 



Ai = hm — log : 

t-oo2t ^^2(0) +^2(0) 



(11) 



is obtained by solving this stochastic differential equation 
by means of a standard method due to van Kampen [^, 
the final analytic expression for Ai reads as 



A,(no,.„,r)4(A-|^ 



where flo = {kR)^,^^ = N {5'^kii) ^j^, 



(12) 



A = I 2ctcV 




— ) + (2fTjV) 



1/3 



and 



2r = 



2^rJo (r^o + o-o) + TTCTo 



(13) 



(14) 



3. ANALYTIC COMPUTATION OF LYAPUNOV 
EXPONENTS 

In the following, we work out time averages of the Ricci 
curvature and of its fluctuations along some analytically 
known unstable periodic orbits of the system described 
by the Hamiltonian 



N N 



i=l 



- qtf 



^(?»+i - qC) 



(15) 



with periodic boundary conditions qiss+x = q\- This sys- 
tem has been introduced by Fermi, Pasta and Ulam in 
their celebrated work on the equipartition properties 
of the dynamics of many non-linearly coupled oscillators. 
Since then, a huge amount of papers have been devoted 
to the study of the link beetwen microscopic dynamical 
properties and macroscopic thermodynamical and statis- 
tical properties of classical many body systems. 

The linear terms in Hamiltonian (|15() can be diago- 
nalized by introducing suitable harmonic normal coordi- 
nates. The latters are obtained by means of a canonical 
linear transformation |lOj| . Denoting the normal coordi- 
nates and momenta by Qk and for A: = 0, . . . , iV — 1 
the transformation is given by 



Qfe(i) 



N 

E 

n=l 



Sknqk{t) , Pk{t) 



N 

E 



SknPk{t) , (16) 



where k — 0, . . . , N ^ 1, and Skn is the ortogonal matrix 
[Toll whose elements are 



1 



sm 



27rfcn 
N 



cos 



27rfcn 
N 



(17) 



n=l,...,N and fc = 0, . . . , TV - 1. The full Hamiltonian 
(|15|l in the new coordinates reads 



HiQ,P) 



2^0 



1 



N-l 



-2^ 

i = l 



Pt+LotQi)+H,{Q) , (18) 



where the anharmonic term is 

N-l 



Hi{Q) = — ^ ujiUjjUJk^^iCijkiQiQjQkQi ■ (19) 



87V 

— 1 

The ijJk = 2sin(7rfc/7V), for fc e {1, . . . , TV - 1}, are the 
normal frequencies for the harmonic case (/i = 0), being 
ujk = tJN-k- By defining 



(— 1)™ for r = mN with m G Z 
otherwise , 



(20) 



the integer-valued coupling coefficients Ciju are explic- 
itly given by 

Cijki = —Ai+j+k+i + Ai+j-k-i + Ai-j+k-i + Ai-j-k+i ■ 

(21) 

By eliminating the motion of the center of mass (which 
corresponds to the zero index), we now easily get the 
equations of motion for the remaining — 1 degrees of 
freedom, which, at the second order, read as 



(3uJr 

2iV 



N-l 



WjUJkUJlCrjklQjQkQl , (22) 



],k,l=l 



for r = 1, . . . ,iV- 1. 

As is shown in Ref.|0|, the equations of motion H22|l 
admit some exact, periodic solutions that can be explic- 
itly expressed in closed analytical form. The simplest 
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ones, consisting of one mode (OM), have only one ex- 
cited mode, which we denote by the index e, and thus 
are characterized by Qj{t) = for j 7^ e. The sohtary 
modes are found by setting Creee — O'^r & {1, . . . , N —1} 
with r ^ e; it is easily verified that this condition is sat- 
isfied for 

_ N N N 2N 3N 

T' Y' T' ~' ~ • ^ 

Thus, for solutions with initial conditions Qj = and 
Qj = for j e, the whole system (I22II reduces to a one 
degree of freedom (and thus integrable) system described 
by the equation of motion 



Qe — ^^eQe ^ 



27V 



(24) 



where Ceeee = 4, 4, 3, 3, 2 for e = iV/4, 3iV/4, N/3, 2N/3, 
N/2, respectively. The harmonic frequencies of the 
modes ^ are uj^ = V2, V2, %/3, V3, 2 for e = N/A, 
3iV/4, N/3, 2N/3, N/2, respectively In order to simphfy 
the notation, in the following, let us set Ce — Ceeee- 
The general solution of (|24(l is a Jacobi elliptic cosine 



,(t) = Acn[rie(t - to) , k] 



(25) 



where the free parameters (modal) amplitude A and time 
origin to are fixed by the initial conditions. The frequency 
fL and the modulus k of Jacobi elliptic cosine function 
depend on A as follows 



rie 



UJeVl + S^ 



SeA^ 



2{1 + SeA^) ' 



(26) 



with Se — f3uj'^Ce/{2N). This kind of solution is periodic, 
and its oscillation period Te depends on the amplitude A, 
since it is given in terms of the complete elliptic integral 
of the first kind K(fc) and in terms of ile by 



Te = 



4K(fc) 

Qe 



(27) 



The modal amplitude A is one-to-one related to the en- 
ergy density e = E/N . In fact, computing the total en- 
ergy (fH^ on the OM solution Qj{t) = SjeQe{t), one finds 



.2/n2 



(28) 



Since at t = the coordinates result (Qe(^o): Pe(^o)) 
{A, 0), by solving the previous equation for A we get 



A 




1/2 



(29) 



This relation allows to express all the parameters of the 
solution (|25|l in terms of the more physically relevant 
parameter e. The period Tg is 



Te 



4K(fc) 



UJe{l + 2peCeV/i 



(30) 



where k = k{e) can be found from l|26() and H29|l . 

In terms of the standard coordinates, the OM solutions 
result 



qn{t) 



1 



N 



dt) 



( 2'Kne\ ( 2'Kne 



(31) 



where e is one of the values listed in H23|l . 

The Ricci curvature along a periodic trajectory, ob- 
tained by substituting Eq. into Eq. jirUjl . is 



kj,{t) = 2^'^-^u:lQl{t) , 
and we can compute its time average kB. as 



6/3 



27« 



(32) 



(33) 



After simple algebra, using standard properties of the 
elliptic functions, we find 



1 

T 



dtQl 



e Jto 



A^ 
Kk^ 



(E+(fc2_i)K) . (34) 



The time averaged Ricci curvature results 
^ „ 12 



Kk^Ce 



1 + 2l3eCe - 1 



[E + (fc2 _ l)K] 



(35) 

where K and E are the complete elliptic integrals of the 
first and second kind respectively, both depending on the 
modulus k which, from (|26|l and (|29|l . is determined by 
the energy density e 



1 



1 + 2(3£Ce 



(36) 



Now, using Eqs. H35I) and (|36|) . and the tabulated val- 
ues for E and K, kji is given as a function of the energy 
density e. In Fig. ^ a comparison is made between kfj 
versus e, worked out for the OM solutions under consid- 
eration, and {kpt) fj^j^ versus e, the average Ricci curvature 
analytically computed in Ref . J9| . 

By definition, the average of the curvature fluctuations 

is 



(37) 



Again, by replacing the microcanonical averages with 
time averages, from Eq. H33|l and after some trivial alge- 
bra, we get 



36/3^w 

iV2 



2, ,4 r 



(38) 



The new term 
Jo 



A^ 



dt cn^(f^et, ^) = / cn'^l^' 



4K 



5 
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FIG. 1: kn versus e, worked out by means of the three sin- 
gle mode solutions identified by the values of e listed in H2.'i|l 
(dotted, dashed and long-dashed lines refer to e = Ai'/4, 3A''/4, 
e = 2Af/3 and e = N jl^ respectively), is compared with 
{kR)fj,^ computed in (continuous line). The agreement is 
very good on a broad range of values of energy density e. 



can be computed by resorting to standard properties of 
the elliptic functions and the result is 



[K(2-5fc2 + 3fc'*) + 2E(2fc2-l)] . (39) 



3Kfc4 

Finally, Eqs. ^ and in JSHJ yield 



192 



(fc2 - 1) + 2(2 - fc 



2^ E n (-E\ 
K 



(1 - 2fc2)2C| 



(40) 



From Eq. H36|) and making use of the tabulated values for 
E and K, equation (|40|l provides the mean fluctuations 
of curvature as a function of e. 

In Fig. 121 a comparison is made between the time aver- 
age of the Ricci curvature fluctuations S'^kp as a function 
of the energy density e, worked out along the OM solution 
that we considered, and (^^fc/})^^ versus e analytically 
computed in Ref. The agreement is very good, thus 
conflrming from a completely new point of view, that 
unstable periodic orbits are special tools for dynamical 
systems analysis; in this case, certain geometric quanti- 
ties of configuration space are surprisingly well sampled 
by UPOs because time averages computed along them 
are very close to microcanonical averages performed on 
the whole energy hypersurfaces. 

Finally, we can compute the Lyapunov exponents as a 
function of the energy density e by inserting Eqs. (|35|l . 



FIG. 2: In this figure we report three curves for Pka versus e 
computed by integrating the curvature fluctuations along the 
three single mode solutions considered in the present paper 
(dotted, dashed and long-dashed lines refer to e = N/A, SN/A, 
e — N/3, 2N/S and e = N/2, respectively), and a comparison 
is made with the same quantity computed in (continuous 
line). Also in this case the agreement is very good. 



(I36|l and (|4()|) into the analytic formulaeElandEl replac- 
ing (fc/j)^£, and {S'^kfi)f^^ by means of the corresponding 
time averages computed above. Fig. ^ shows that the 
overall agreement between our analytic results, the ana- 
lytic results from ^ and the results obtained by numer- 
ical integration of the tangent dynamics, is very good. 
The agreement is globally very good because at high en- 
ergy density our results are really very close to the other 
mentioned ones, and at low energy density the discrep- 
ancy does not exceed - at worst - a factor of 2 on a range 
of many decades of energy density and with the use of 
only one unstable periodic orbit! 



4. CONCLUDING REMARKS 

In conclusion, we have found that some global cur- 
vature properties of the configuration space manifold 
- whose geodesies coincide with the trajectories of an 
Hamiltonian system - are well sampled by unstable pe- 
riodic orbits. Then, since the averages of these curva- 
ture quantities enter an analytic formula to compute the 
largest Lyapunov exponent, unstable periodic orbits can 
be used also to compute Lyapunov exponents through 
the time averages of the same geometric quantities. In 
the present work, this result has been obtained in the 
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case of the FPU-/3 model for which the analytic expres- 
sion of some unstable periodic solutions of the equations 
of motion are known. The outcome of this computa- 
tions is in very good agreement with those reported in 
Ref.0| on Ai for the same model. Of course, it would be 
very interesting to perform similar computations also for 
other models. Finally, from a new point of view, that of 
the Riemannian geometric theory of Hamitonian chaos, 
we confirm that unstable periodic orbits seem to have 
a special relevance among all the possible phase space 
trajectories of a nonlinear Hamiltonian system. 



10"^ 10"^ 10"' 10° 10' 10' 10^ 10' 10° 



FIG. 3: This figure shows the largest Lyapunov exponent Ai 
obtained by integrating the suitable geometric quantity along 
the three single mode solutions considered in the present pa- 
per, plotted vs. e. Dotted, dashed and long-dashed lines refer 
to e = iV/4, 3iV/4, e = N/Z, and e = N/2, respectively. 
Continuous line refers to the Lyapunov exponent computed 
in 0. The full circles are the values for Ai computed by nu- 
merical integration. The agreement is again very good on a 
broad range of e values. 
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